Real-time dynamics at finite temperature by DMRG: 
A path-integral approach 
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We propose a path-integral variant of the DMRG method to calculate real-time correlation func- 
tions at arbitrary finite temperatures. To illustrate the method we study the longitudinal autocor- 
relation function of the XXZ-chain. By comparison with exact results at the free fermion point 
we show that our method yields accurate results up to a limiting time which is determined by the 
spectrum of the reduced density matrix. 
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The density-matrix renormalization group (DMRG) 1 
is today a well established numerical method to study 
ground state properties of one-dimensional quantum sys- 
tems. Within the last few years the DMRG method 
has been generalized to allow also for the calculation of 
spectral functions and, quite recently, to incorporate di- 
rectly real-time evolution^ The most powerful variant 
of DMRG to calculate thermodynamic properties is the 
density-matrix renormalization group applied to transfer 
matrices (TMRG). This method has been proposed by 
Bursill et al. A and has later been improved considerably. 5 
The main idea of TMRG is to express the partition func- 
tion Z of a one-dimensional quantum model by that of an 
equivalent two-dimensional classical model obtained by 
a Trotter-Suzuki decomposition. 6 Thermodynamic quan- 
tities can then be calculated by considering a suitable 
transfer matrix T for the classical model. The main ad- 
vantage of this method is that the thermodynamic limit 
(chain length L — > oo) can be performed exactly and 
that the free energy in the thermodynamic limit is deter- 
mined solely by the largest eigenvalue of T. The TMRG 
has been applied to calculate static thermodynamic prop- 
erties for a variety of one-dimensional systems including 
spin chains, the Kondo lattice model, the t — J chain and 
ladder, and also spin-orbital models£2*2i 

The Trotter-Suzuki decomposition of a one- 
dimensional quantum system yields a two-dimensional 
classical model with one axis corresponding to imaginary 
time (inverse temperature). It is therefore straightfor- 
ward to calculate imaginary-time correlation functions 
(CFs) using the TMRG algorithm. Although the results 
for the imaginary-time CFs obtained by TMRG are very 
accurate, the results for real-times (real-frequencies) 
involve large errors because the analytical continuation 
poses an ill-conditioned problem. In practice it has 
turned out that the maximum entropy method is the 
most efficient and reliable way to obtain spectral func- 
tions from TMRG data. The combination of TMRG and 
maximum entropy has been used to calculate spectral 
functions for the AAZ-chain 1Q and the Kondo-lattice 



modelfSi However, this method involves intrinsic errors 
due to the analytical continuation which cannot be 
resolved. 

Here we propose a method to calculate directly real- 
time CFs at finite temperature by a modified TMRG 
algorithm thus avoiding an analytical continuation. We 
start by considering the two-point CF for an operator 
O r (t) at site r and time t 
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where j3 is the inverse temperature. Here we have used 
the cyclic invariance of the trace and have written the 
denominator in analogy to the numerator. For a Hamil- 
tonian H with nearest-neighbor interactions we can use 
the Trotter-Suzuki decomposition 
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where e = (3/M so that the partition function Z = 
exp(— (3H) becomes 
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With it — > t in Eq. Q and inserting the identity operator 
at each imaginary time step one obtains directly a lattice 
path-integral representation for the imaginary time CF 
(O r (r)O (0))Altt 

The crucial step in our new approach for real times 
is to introduce a second Trotter-Suzuki decomposition of 
exp(— iSH) as in Eq. @ with S — t/N . We can then 
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define a column-to-column transfer matrix 

T~2N,M = (n, 2^3,4 •• ■7"2Ar-l,2Af)(T2,37"4,5 ••■ T2Af,2M+l) 
(v2M+l,2M+2 ■ ■ ■ V2M+2N -1,2M+2n) 
(v2M+2,2M+3 ■ ■ ■ V2M+2N,2M+2N+l) 
(v2M+2N+l,2M+2N+2 1 1 ■ V2M+4jV-l,2M+47v) 
(^2M+2iV+2,2M+2iV+3 1 1 ' W2M+47V,l) (4) 

where the local transfer matrices have matrix elements 
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and v is the complex conjugate. Here i = 1, • • • , L is the 
lattice site, fc = 1, • • • , 2M (I = 1, ■ ■ ■ , 2 A) the index of 
the imaginary time (real time) slices and Sw^ denotes a 
local basis. The denominator in Eq. can then be rep- 

L 12 

resented by Tr(X,^ M ) where N,M,L — > oo. A similar 
path-integral representation holds for the numerator in 
(JTJ. Here we have to introduce an additional modified 
transfer matrix T2n,m{0) which contains the operator O 
at the appropriate position. For r > 1 we find 
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Here [r/2] denotes the first integer smaller than or equal 
to r/2 and we have set T = T^jv.m- In the second step 
we have used the fact that the spectrum of the column- 
to-column transfer matrix T is gapped for j3 < oo. 5 
Therefore the trace is reduced to an expectation value 



where 1^^) 



(vPq | are the right- and left-eigenvectors of 



the non-hermitian transfer matrix T belonging to the 
largest eigenvalue Aq. A graphical representation of the 
transfer matrices appearing in the numerator of Eq. © 
is shown in Fig. ^ 

It is worth to note that Eq. © has the same struc- 
ture as the equation for the calculation of imaginary time 
CFs&iS Only the transfer matrices involved are differ- 
ent. For practical DMRG calculations the parameters 
e, S are fixed and the temperature (time) is decreased (in- 
creased) by increasing M (N). This is achieved by split- 
ting T, T(0) into a system and an environment block (see 
Fig. . Adding an additional r-plaquette at the lower 
end of the system block then decreases the temperature 
whereas adding a w-plaquette at the upper end leads to an 
increase in time t. The extended system block is renor- 
malized by projecting onto the Z largest eigenstates of 
the reduced density matrix ps = Tie\^q)(^o\- Starting 
from an existing TMRG program our real-time algorithm 
requires only the following minor changes: In addition to 
the T-plaquettes also w-plaquettes are present. Therefore 
the transfer matrix and the density matrix become com- 
plex quantities so that complex diagonalization routines 
are required. 
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FIG. 1: Transfer matrices appearing in the numerator of 
Eq. © for r > 1 with r even. The 2 black dots denote the 
operator O. Each transfer matrix consists of three parts: A 
part representing exp(—(3H) (vertically striped plaquettes), 
another for exp(itH) (stripes from lower left to upper right) 
and a third part describing exp(— itH) (upper left to lower 
right). The transfer matrices T,T(0) are split into a system 
(S) and an environment block (E). 



In the remainder of this letter we want to study, as 
a highly non-trivial example, the longitudinal autocorre- 
lation C(t,T) = (S${t)S$(0)} at temperature T for the 
XX Z chain with Hamiltonian 
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where S = 1/2, J > and A > 0. Although the model is 
integrable, exact results for C(t,T) are available only at 
the free fermion point A = 0. Quite interestingly, C(f , T) 
shows a non-trivial behavior even at infinite temperature 
due to the quantum nature of the problem^i 

For the autocorrelation both S z operators are situated 
in the same transfer matrix T(S Z ,S Z ) so that Eq. © 
reduces to 
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It is important to note that once the blocks necessary 
to construct T ,T(S Z , S z ) at a given time t are known, 
the correlation function for arbitrary distance r with 
t fixed can be simply calculated by additional matrix- 
vector multiplications (see Eq. ©). We start with the 
case T = oo where T ,T(S Z , S z ) do not contain any r- 
plaquettes. This limit directly addresses the essential 
feature in our approach, namely the transfer-matrix rep- 
resentation of the time evolution operator. In Fig. [21 re- 
sults for A = and A = 1 are shown where the number 
of states kept in the DMRG varies between Z = 50 and 
Z = 400. In the case A = the autocorrelation func- 
tion can be calculated exactly by mapping the system to 
a free spinless fermion model using the Jordan- Wigner 
transformation. For arbitrary temperature T the result 
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FIG. 2: Longitudinal autocorrelation for A = and A = 1 
(inset) at T = oo where the number of states kept varies 
between Z = 50 and 400 and 5 — 0.1. The exact result is 
shown for comparison in the case A = 0. 

can be expressed as C(t, T) = {D x (t) + D 2 (t, T)} 2 with 
„ , . J (Jt) , s i f 1 sin jet , e ,„ 

(9) 

where Jo is the Bessel function of order zero>i£ At T = 
oo this reduces to C(t, oo) = Jg(Ji)/4. Our numerical 
results agree with the exact one up to a maximum time 
which is determined by the number of states kept in the 
DMRG. For the maximum number of states, Z = 400, 
considered here we are able to reproduce the exact result 
up to Jt ~ 12. For A = 1 no exact result is available to 
compare with, however, the A = case suggests that the 
results are trustworthy at least as long as they agree with 
the data where a much smaller number of states is kept. 
The numerical data with Z = 400 should therefore be 
correct at least up to Jt ~ 10. We also checked the A = 1 
data for small t by comparing with exact diagonalization. 
For TMRG calculations the free fermion point is in no 
way special and the algorithm is expected to show the 
same behavior at this point as for general AiS In the 
following we will therefore concentrate on A = where 
a direct comparison with exact results is possible. 

For small t the error in the numerics is entirely domi- 
nated by the finite Trotter-Suzuki decomposition param- 
eter 8. It is therefore possible to enhance the accuracy by 
decreasing 8 as shown in Fig.JSJa). As expected, the error 
is quadratic in 8 but fortunately with a rather small prc- 
factor ~ 10 -3 . For small 8 more RG steps are necessary 
to reach the same t. Interestingly, the breakdown of the 
algorithm does not depend on the number of RG steps. 
For different 8 but fixed Z it always occurs at about the 
same time t c . Furthermore the breakdown is always a 
very rapid one, i.e., for times considerably larger than 
t c the errors become arbitrarily large. This suggests that 
there is an intrinsic maximum time scale set by the prob- 
lem itself. This is supported by Fig. EJb) showing a rapid 



increase in the number of states Z n necessary to keep the 
error below 10~ 3 . To understand this behavior we have 
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FIG. 3: (a) A e = \C(t, oo) e xact — C(t, oo)tmrg\ at t = 4 as a 
function of S with Z — 100 states kept. The line corresponds 
to a fit with A e = 3.05 • 10~ 3 <5 2 . (b) Number of states Z n 
required to keep the error below 10 -3 as function of t for 
T/J = oo,10,l,0.1. (c) Largest 100 eigenvalues Ai of ps at 
T — oo calculated exactly. The inset shows the discarded 
weight 1 - Si=i A,:- 

calculated the spectrum of the reduced density matrix ps 
exactly for T — oo and A = 0. The result is shown in 
Fig- Etc)- For small t the spectrum is decaying rapidly 
so that indeed a few states are sufficient to represent the 
transfer matrix T accurately. At larger time scales, how- 
ever, the spectrum becomes dense. This means that the 
number of states needed for an accurate representation 
starts to increase exponentially in agreement with our 
numerical findings shown in Fig|3Jb). The breakdown 
approximately occurs when the discarded weight defined 
by 1 — J2i=i A»> where Aj are the Z largest eigenvalues of 
ps, becomes larger than 10 -3 (see inset of Fig.^c)). The 
long-time asymptotics of the autocorrelation function at 
T = oo is therefore not accessible within our method. 

Next, we consider finite temperatures < T < oo. 
Results for A = and T/J = 10,1,0.1 are shown in 
Fig. 0] According to Eq. 10, C(t,T) now acquires also 
an imaginary part which is shown in the insets of Fig. 0] 
For T I J — 10 and T/J — 1 the results look qualitatively 
similar to the T = oo case. As shown in Fig. I^b) the 
number of states needed to obtain the same accuracy 
at a given time increases with decreasing temperature. 
This is easy to understand because the Hilbert space for 
T increases when adding additional r-plaquettes. The 



4 



0.25 
- 0.2 

r 

i 0.15 

) 

: o.i 

0.05 




T/J = 


10 














(<(0) 










- 


- \ 














- \ 


V 










200 J 

exact 








1 


i 


100 ^ 
1 , 1 


~— ^ 300 / '\^ /- 

- y 400 --^r^ - 
i.i, 




\ o 


2 


4 


6 8 


10 12 












Jt 

100^_ 


300 >— 
400 




1 




l 


l 


50—5^^ 


exact 

200 ^ ^ s 



8x10 
6x10 
4x10 
2x10 


-2x10 



-3 




0.5 
0.4 
0.3 
0.2 



T/J = 0.1 




















o 










N w 










N v 

A ¥ 




\ 50 X/ 


200 


exact 






\ / 100 — 


1 


40</ 300^ . 
1,1,1, 


- \ 


2 4 6 8 


10 


12 14 16 






~\ 50 X^\ • • • 




400 » 






100 ^ 




300"^ exact 






1,1,1 


1 


1,1, 



0.1 



10 



FIG. 4: C(t,T) for A = at T/J = 10,1,0.1. The main 
figures show the real, the insets the imaginary parts. 



breakdown of the algorithm also looks similar to the T = 
oo case and we again find an exponential increase in Z n 
at larger times. As we have chosen e = 5, the number 
of r-plaquettes M is much smaller than the number of 
v-plaquettes N for t ~ 10 where the calculations with 
Z = 300,400 states fail. The spectrum of the reduced 
density matrix ps at these time scales will therefore look 
very similar to the one for T — oo shown in Fig. IHIc). 
We therefore conclude that it is again the spectrum of 
ps which sets the limiting time for our calculations. 

For T/J = 0.1, however, we find a different behav- 
ior. Instead of a rapid breakdown with arbitrarily large 
deviations for t > t c we find that the results for all Z 
remain relatively close to the exact one over the entire 



time scale investigated (the data for Z = 50, 100 in Fig.0] 
are depicted only up to times marking the start of de- 
viations). We also see from Fig. Ofb) that the func- 
tional form of the increase in Z n is now different from 
the cases T/J = oo, 10, 1. Whereas it becomes exponen- 
tial in the latter, it is more or less linear for T/J = 0.1 
from Z n = 200 to Z n = 400. This can be understood 
as follows: For a transfer matrix T consisting only of 
r-plaquettes, the spectrum of ps is exponentially decay- 
ing. This characteristic should be preserved as long as 
the number of v-plaquettes is not much larger than the 
number of r-plaquettes. A fundamental failure of our 
approach due to a dense spectrum of ps should only oc- 
cur for N ^> M. By considerably increasing Z it should 
therefore be possible to access much larger time scales at 
low temperatures in future large scale numerical studies. 

To conclude, we have presented a numerical method to 
calculate real-time correlations in one-dimensional quan- 
tum systems at finite temperature. The method is based 
on a DMRG algorithm applied to transfer matrices. As 
essentially new ingredient it involves a second Trotter- 
Suzuki decomposition for the time evolution operator. 
To test our approach we have calculated the autocorre- 
lation function for the XXZ-chain both at infinite and 
finite temperatures. For T — oo we have established that 
reliable results can be obtained up to a maximum time 
scale t c where the spectrum of the reduced density ma- 
trix ps becomes dense. For low T we have shown that 
the algorithm does not show a rapid breakdown contrary 
to the high-T case and have argued that the fundamen- 
tal problem of ps becoming dense is less severe. A huge 
advantage of our approach compared to other methods is 
that once the blocks necessary to construct T, T{0) at a 
given time t are known, the CF can be evaluated at time 
t for arbitrary distances r between the operators O by 
simple matrix-vector multiplications. As our approach is 
working directly in the thermodynamic limit it is possible 
to obtain highly accurate results even for large distances 
as has been demonstrated for the static case in Ref. 
This will be exploited in a forthcoming publication^ for 
a detailed study of time (t < t c ) and space dependent 
CFs in the XX Z chain. 
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